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CJH 1 Abstract. 

^N , A numerical simulation of vibrational excitation of molecules was devised, and 

f-*. \ used to excite computational models of common molecules into a prescribed, pure, 

normal vibration mode in the ground electronic state, with varying, controlable energy 

content. The redistribution of this energy (either non-chaotic or irreversible IVR) 
qq ' within the isolated, free molecule is then followed in time with a view to determining the 

^»0 • coupling strength between modes. This work was triggered by the need to predict the 

general characters of the infrared spectra to be expected from molecules in interstellar 
^^ ■ space, after being excited by photon absorption or reaction with a radical. It is found 

VO . that IVR from a pure normal mode is very "restricted" indeed at energy contents 

of one mode quantum or so. However, as this is increased, or when the excitation 
rSZi • is localized, our approach allows us to isolate, describe and quantify a number of 

interesting phenomena, known to chemists and in non-linear mechanics, but difficult 
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to demonstrate experimentally: frequency dragging, mode locking or quenching or, 
still, instability near a potential surface crossing, the first step to generalized chaos as 
— j \ the energy content per mode is increased. 
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1. Introduction 

The work to be described here is an effort to simulate the excitation of pure vibrational 
modes in a single molecule, then follow their evolution in time, or relaxation paths, with 
a view to estimating coupling strengths between modes and determining bottlenecks in 
the way to the "final" molecular vibration state, that is, before possible subsequent 
infrared emission, which occurs much later in general. 

This effort was triggered precisely by the need to predict, and hence be able to 
interpret, the observed IR (infrared) emission from isolated interstellar molecules which 
have been excited by radiation absorption or recombination of radicals (see Williams 
and Leone 1995J). Such excitation usually leads to electronic state transitions, internal 
conversion and intersystem crossing, and, at least, mixing of vibrational states. This 
was masterly described in the pedagogical overview written by Nesbitt and Field |199fij . 
to which the present work owes much. 

Already in the 80's, Zewail and collaborators |1984aj . 1984b had illustrated 
experimentally, and analyzed theoretically irreversible IVR, involving a quasi-continuum 
of interacting states, but also pointed out the importance for chemistry of non-collisional, 
non-chaotic, multilevel energy flow, involving only a few coupled states, as evidenced by 
quantum beats observed in the fluorescence of the excited molecule. They also pointed 
out the analogy with systems of coupled mechanical oscillators. 

At about the same time, it was realized that relaxation in an isolated molecule (non- 
collisional relaxation) did not always quickly end up in "thermodynamic equilibrium", 
as if the molecule were embedded in a thermal bath: rather, it was necessary to 
acknowledge the existence of mode selectivity (see Sewell et al [1992 ), which may forbid 
certain regions of molecular phase space and lead to "restricted" and quasi-periodic 
flows, as opposed to "facile" , chaotic or downright ergodic motion, which characterizes 
thermal equilibrium (see Noid et al [1981J; Brumer [19^5)- Chemists were relieved to 
notice that this need not contradict the celebrated RRKM treatment of uni-molecular 
reactions: the latter is most generally valid since it does not in effect require complete 
ergodicity (Sewell at al. |1992j ). and since, at sufficiently high excess energy density in the 
molecule and depending on the extent of overlapping nonlinear resonance interactions 
(Sibert et al [1984J, Ford[1992 ), ergodicity does set in, as assumed in the first place. 

Because of its potential use in selective chemistry, the transition between quasi- 
periodic and chaotic regimes has since been the subject of considerable interest and 
numerous experiments (e.g. Sibert et al [1984J, Sewell et al [1992J, Deak et al [1999J, 
Stromberg et al |2002]). The experimental approach to these phenomena requires 
sophisticated apparata, including picosecond spectroscopy and cold molecular beams. 
Also, for experimental convenience, excited electronic states are privileged over their 
ground state. Unfortunately, the former lack the desired IR information. Finally, 
the information about the excited states is mostly deduced indirectly from whatever 
diagnostic means are available. 

Here, we are only interested in the quasi-periodic, or coherent regime, in a limited 
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range of molecules (medium-size hydrocarbons) and take advantage of the recent 
development of highly performing chemical computational software to simulate the types 
of molecules we are interested in and their behaviour upon specified excitation into 
higher vibration levels in the ground electronic state, far from the dissociation limit. 

Moreover, the molecular behaviour we are interested in is strongly affected by 
anharmonocity, and so depends heavily on the energy content of the target excited 
state. Again, the latter parameter is much easier to control computationaly than 
experimentally, as shown below. 

Finally, the molecules of astrophysical interest may be composed of as many as a 
few hundred atoms, and so are very difficult to produce and isolate for excitation and 
detection, making inevitable the resort to numerical modeling. 

For all the imperfections of chemical codes, they are quite convenient (as already 
noticed long ago: see Noid et al |1981| ). because they provide the possibility of collecting 
directly the desired information with very high time resolution (10 _1 fs): in each 
experiment, all molecular properties are collected as as a function of time, including all 
bond lengths and orientations, as well as global properties such as kinetic and potential 
energies, electric dipole moment, etc. 

A key step in carrying out our project was the design of a computational procedure 
to energize a single normal mode in a chemical model of the real molecular structure, 
so as to be able to isolate and define clearly individual relaxation paths. This produced 
sparse spectra and simple beating patterns, from which to deduce coupling energies. 

The simulation of relaxation from high energy levels also uncovered chaotic 
behaviours, difficult to observe experimentally, but well known in non-linear mechanics, 
when exploring the phase space of non-integrable systems, starting with bifurcations 
to end up in various forms of chaos. Here are demonstrated, for instance, multiple 
resonances, localized modes, frequency drag, mode locking (synchronization) and mode 
quenching, some of which are seldom treated in the current litterature. 

In the following, Sec. 2 is a general reminder of the behaviour of coupled linear 
oscillators, as an introduction to the formalism, together with the non-linear corrections 
which were needed to explain the observations described in this work, on common 
molecules, even at moderate energy contents (Av = 1 or 2). Of course, all of this is 
best done quantum mechanically (Nesbitt and Fielcl[1996j, Gruebele and Wolynes[2004j, 
Wu 2005J), but the parallel is so perfect that the classical treatment, including Classical 
Dynamics for the motion of nuclei (Born-Oppenheimer approximation), can deliver much 
of the right physical picture (e.g. Noid et al [1981J, Kuharski et al [1988 ) (except for 
quantization and zero-level energy which do not matter here) with only a very light 
mathematical apparatus (Herzberg.ljMSJ. 

Section 3 describes a new procedure for the excitation of single modes and the 
effects of various parameters defining this procedure. Sections 4 to 9 apply the above to 
a simple computer-simulated molecule, C 2 H 4 . Several excitation and relaxation routes of 
interest are considered. Because of the relative sparseness of its vibrational spectrum, 
it is possible to isolate interesting molecular phenomena and determine some of the 
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coupling factors involved. Section 10 is an extension to a much larger molecule (101 
atoms) with a correlatively much denser spectrum. 

A general conclusion drawn from these observations is the weakness of most mode 
couplings and the unequal partition of the excitation energy among the modes at the 
end of the relaxation path. Only at higher energy contents, the higher the larger the 
molecule, are chaos and better equipartition, observed (Papoular 20U2]). but that is not 
considered here. 

2. Coupled non-linear oscillators 

In the ground vibrational state the potentials which determine the nuclear motions are 
parabolic in first approximation. As a consequence, the motion of each atom can be 
broken down into 3N-6 independent, collective , normal modes, where N is the number 
of nuclei in the structure. 

However, as soon as the molecule is set in motion, deviations from the quadratic 
dependence of potentials on excursions of nuclei from their equilibrium positions 
(anharmonicity) changes the normal mode frequencies and couples all these oscillators 
together, so that energy can flow from one to the other, more or less easily, depending 
on their coupling coefficients and frquency separations. The ways they do so are the 
subject of this study. 

Consider two oscillators with frequencies V\ and u 2 , coupled by a perturbation term, 
v, in the equilibrium potential. In quantum mechanics, this is an off-diagonal term in 
the Hamiltonian matrix; in electric circuit theory, it is expressed by the impedance of an 
electric component common to both coupled circuits. In mechanics, it is the stiffness, k, 
of a spring linking two pendula; for instance, if m is their common mass, the equations 
of motion are: 

£i + wfxi + ^(xi-ar 2 ) = 

X 2 + Lo\x 2 ~ £{X! ~ X 2 ) = 0, 

where uo = 2-kv. In the important particular case of resonance, v\ = v 2 = vq, these 
equations are easily separated into two independent harmonic motions by transforming 
variables into X a = x\ + x 2 and X b = x\ — x 2 , giving 

X a + uj 2 X a = 

x b + K 2 + fjx b = o 

Obviously, the natural frequencies of these so-called normal modes are 

v a = v ;v b = v (l + ^\v = ^- o . (1) 

One of the new frequencies coincides with the common, unperturbed one, while the 
other is slightly blue shifted. This corresponds to the quantum classical repulsion of 
crossing (resonant) energy levels. For sufficiently weak coupling, v b ~ vq(1 + — ). These 
equations perfectly illustrate, e.g., the symmetric and anti-symmetric scissor modes of 
C2H4, which is treated below: x\ )2 representing the two HCH angles, while X ajb are the 
two normal modes. 
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Returning now to the general, out-of-resonance, case, the solutions are of the form 
X\ (t) = Acosu) a t + BcOSLUbt 
X2(t) = Ccosu a t + Dcosubt, 
where 

v a , h = {^ + vfavtfl* T [(^) 2 + v^v^yl* (2) 

and A, B, C and D are constants. Each oscillator is a linear combination of modes 
a, b, with relative amplitudes depending on initial conditions. 

The most telling case, then, is when oscillator 1 is initially excited alone. A simple 
trigonometric transformation of the equations for X\p reveals that, when the coupling 
is not too strong, the motion of each oscillator can also be viewed as sinusoidal, with 
common frequency 

«a±ja+* (3) 

and amplitude-modulated with frequency 

fi=[I^)! + ^]V2. (4) 

This modulation implies that a corresponding amount of energy is periodically 
exchanged between the oscillators, as illustrated in fig. 1 (cf. Zewail'z quantum beats). 
It turns out that the average fraction of energy that is delivered to oscillator 2 is 

e=[l + ( !a ? a ) 2 ]- 1 - (5) 

Equation (5) predicts that the maximum amount of energy that can be transferred 
from 1 to 2 increases as the coupling increases and the frequency separation decreases, 
to reach unity when their ratio becomes infinite (i.e. at resonance). Also, the transfer 
frequency (4) increases with both parameters. Thus, even in steady state, there is no 
equipartition in general. In a quantum mechanical context, the squared amplitude of 
a mode i should be viewed as the probility of finding the system in state % at time t; 
complete transition from state 1 to state 2 (e = 1) can only occur in case of resonance 
and then only periodically, as in Rabi's oscillations (see Cohen- Tannoudji et al |1987j ). 
When the coupling becomes too strong, however, the separation into two modes no 
longer makes sense. 

Now, a molecule can be considered as a system of oscillators in different ways. 
In ethylene, for instance, the two opposite pairs of C-H bonds can be considered 
as resonant angular oscillators coupled by the C-C bond, as we do in Sec. 4. In 
a different analogy, when the energy content drives the system out of the parabolic 
region of the potentials, so the various normal modes determined at low energies can 
no longer be said to be independent, they can instead be considered as oscillators 
linked by perturbative couplings. One can then determine the order and extent of 
such perturbations as a function of energy content. Here, however, the coupling factor 
also depends on the relative orientations of the vectors describing the nuclear motion in 
each mode (symmetry selection rules), a constraint which must be taken into account 
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Figure 1. Envelopes of the peak amplitudes of excursions X\ t 2 of two coupled 
oscillators or normal modes, illustrating (4) and (5). Initially mode 1 alone is excited. 
Both amplitude and periodicity, T — 2n/D,, of the curves carry information on both 
the frequency separation and extent of coupling between the two oscillators. 



when considering the relaxation ladder and considerably restricts IVR (Felker and 
Zewail [1984b]). For instance, in the case of a plane molecule such as benzene or ethylene 
(studied below), out-of-plane CH bends are extremely difficult to excite through the 
"planar" modes. 

Things are further complicated by the anharmonic dependence of couplings and 
mode frequencies upon energy content. However, this also gives rise to interesting new 
phenomena, such as frequency dragging and mode locking, observed below. In order to 
understand the occurrence of such phenomena, the equations above must be modified 
to account for deformations of the potential surfaces. Our approach is: 

- anharmonicity of nth order is approximately simulated by multiplying io- h by 
(1 — ax™), where, as before, x is an oscillation amplitude while n defines the order of 
anharmonicity; the dimension of a is that of x~ n so as to make ax n dimensionless; this 
implies a dissociation limit at Xd = a 1 /™; 

-resonances of second order (u± = i/ 2 + z/ 3 , including Fermi resonance as a special 
case) are approximately simulated by adding to each equation a term b(xi — x 2 ) 2 . 

At that point, of course, the equations must be solved numerically for particular 
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values of the parameters. This was done at every step of this work, so as to understand 
the physical origin of features exhibited by real molecules and described below. The 
individual role of each non-linear term is highlighted although they both stem from the 
same deformation of potential surfaces with increasing energy content and, therefore, 
are not separable. 

3. Energizing a normal mode 

A particular case of the above is the feeding of energy from an external source (e.g. a 
single-mode laser) into a single oscillator (a molecular normal mode) of frequency Uq 
and initially at rest. Classically, it is well known that, in response to an an applied 
resonant perturbation, the oscillator amplitude increases as t and its energy as t 2 , which 
corresponds to Q = in (4). In order to excite one mode of a model molecule in this 
way, it is in principle necessary to devise a means to perturb all bonds in tune, at the 
right frequency and according to the given mode geometry. This is presently almost out 
of reach in most practical cases because of software limitations. However, we show here 
that the desired result can be approached by periodically applying short perturbative 
pulses to a few specific bonds at the same time. If the pulse period is T = — , this is 
equivalent to a harmonic series of frequencies ^ (n integer > 0). It turns out that, if u is 
slightly smaller than ujq, the procedure nicely takes care of the progressive, anharmonic, 
decrease of Uq as energy is fed into it, and that harmonics are not significant. 

Let us now quantify these considerations to optimize the result. This is done with 
the help of fig. 2, where time is in abscissa and bond length (or oscillation amplitude) 
in ordinate. At times t = nT, one or more particular bonds are set to the same (small) 
value, qo, of excursion from equilibrium in length or orientation. Let ujqT = 2ir(l — 5). 
In between two successive pulses at t — and t = T, e.g., the free "motion" of a bond 
is 

Q\ = qicos(uj t + 0i ). 

After the 2nd pulse, the motion is similarly described by 

Ql = q 2 cos(uj T + 02 ) 

with r = t — T. The procedure imposes 

Qi(0)=Q 2 (0) = g , 

so there is, at t — T, a discontinuity between Q\{T) = gicos0i and (^2(0) = q^cos^- 
However, their derivatives must be continuous, i.e. 
qisin((fii + 2tt(1 — 5)) = q 2 sin(f)2. 
Combining these equations gives 

sin[<j>i + 2vr(l - 5)} sin(<f>i - 2ir5) 
tan(p2 = = . (6) 

COS01 COS(pi 

Clearly, phase 0» changes all along, tending to a limit 

0^ = ^(5-1/2), for 5 >0, (7) 

0; = 7^+1/2), for <5 <0. (8) 
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Figure 2. One period of the energizing process, illustrating Sec. 3. Q is the free 
motion of the oscillator that is being energized, between two succeeding, externally 
imposed, excursions, go- 



Correlatively, amplitude q(n) tends in time towards a corresponding upper limit 
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depending on whether S > or< respectively. In fig. 2, this asymptotic limit 
(saturation) corresponds to a situation where, because of the increase in amplitude, 
both ends of the arch of sinusoid have ordinates go- It a l so appears that qi — > oo as 
5^0, but this has a price in that each pulse then delivers a correspondingly smaller 
amount of energy to the target mode, thus slowing the process. As | S | increases, the 
saturation limit is reached earlier but is less energetic. 

As the energy to be deposited may reach a few tens of eV (IR photon energies), 
anharmonicity cannot be ignored: during the energizing process, ujq decreases and S 
increases. By way of an illustration, put 



8( n ) = 5 + aE™, E n = q 2 n /2. 



(10) 



Here, S is the dimensionless slip in frequency, n is the number of elapsed periods, 
m is the order of anharmonicity, E n is proportional to the oscillator energy at time 
t = nT and the dimension of a is such as to make aE™ dimensionless (this convention is 
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Figure 3. Computed examples of energizing curves. Symmetric initial detunings 
between the source and receiver oscillators yield different limit energies: negative 
detunings are favored because of anharmonic red-shift. As defined by eq. (10), E 
is proportional to the oscillator energy at time t. 



slightly different from that of Sec. 2, but carries the same information as to the degree 
of anharmonicity) . 

Figure 3 shows the evolutions for two symmetric initial detunings, 5q = ±0.05, and 
the same non-zero anharmonicity, a = 0.5. Negative initial detunings ultimately deliver 
more energy to the target mode: this is because resonance is approached more closely as 
the effective detuning decreases in the process; in the limit when the mode synchronizes 
with the perturbation, its energy tends to 
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if m is constant. In the absence of anharmonicity, of course, both curves coincide, 
and both limiting deposited energies decrease as the initial detuning increases. The 
optimum value of the latter is highly sensitive to a and m, and must be determined 
by cut-and-try; but a few trials show that, in general, a few percent negative initial 
detuning {ujqT > 2n) yield enough energy after a few tens of energizing cycles. Thus, in 
the case of fig. 3, 5q = —0.05 yields a limiting amplitude gain — = 8; this corresponds 
to an energy gain of 64 w.r.t. the potential energy perturbation induced by the initial 
perturbation, q , of the equilibrium bonds. 
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As to this initial perturbation, it should match as closely as possible the geometry 
of the target mode. While this may be easy for small molecules (see Sec. 4 ff), it is 
hopeless for large ones (Sec. 9). In that case, however, one can select a few bonds 
of the adequate type, and perturb them in the right way ( stretch, in-plane or out-of- 
plane). Thus, it was found in this study that an initial C-C stretching perturbation, for 
instance, is quite helpful in exciting many skeleton and functional modes, albeit with 
reduced efficiency 

When the periodic perturbation does not perfectly match the desired mode 
geometry, what is excited is a localized mode (a "bright" state), i.e. a linear 
superposition of normal modes, which may blur the relaxation path. However, the 
frequency selectivity of the energizing procedure is often sufficient to eliminate residual 
undesirable modes. Note that the situation is quite similar to that of the excitation 
of modes in a microwave cavity by application of a voltage to a suitably shaped and 
oriented metallic probe protruding into the cavity. 

Following the remark at the beginning of this section, the selectivity of the proposed 
procedure can be improved by increasing the number of applied perturbations (g ) per 
cycle, for instance including —go at t — T/2, which eliminates all even harmonics. 

Once an energizing curve has been "experimentally" obtained, it is possible to 
determine the number of cycles necessary to deposit exactly the desired number of 
quanta in the target mode for further study; for instance, to determine the anharmonicity 
parameters, or mode couplings as a function of energy content. 

Note again that, between two succeeding applications of go, the system oscillates 
freely; the chemical numerical simulation codes used here to implement the whole scheme 
are Classical Dynamics codes associated with either the semi-empirical code AMI for 
small molecules (using quantum mechanics to determine electronic configurations at 
each step of the dynamics), or the Molecular Mechanics code MM2 (using experimentally 
adjusted force fields) for large molecules, as embedded in the Hyperchem 5.0 package 
(Papoular2002j). An excellent overview of these codes, as well as references to the 
original papers can be found in Hinchliffc 2003J. 

In the following computations, the ambient temperature is and all dynamics are 
computed by steps of 1 fs. 

4. Case study: C 2 H 4 , coupling energy 

C2H4 has 12 normal modes between 800 and 3200 cm -1 . Two of them, the symmetric 
(sym) and anti-symmetric (a-s) scissor modes, provide a good example of the different 
ways coupled-oscillators theory can be applied. In the first, both CH bond angles 
oscillate in phase, and in the second, with opposite phases. These modes can be 
considered as two oscillators, isolated at low energy contents, which are coupled by 
anharmonicity at higher energy levels. However, the symmetry of the system also 
allows us to view the two pairs of CH bonds as two identical oscillators, beeing coupled 
physically by the CC bond. As a result, the motion of the system will break down into 
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2 modes, one at the initial frequency of the isolated oscillators, the other slightly blue- 
shifted. Normal mode analysis places them at 1388 and 1412.5 cm" 1 , to the accuracy 
of the present simulation. Using (1) with v a — vq = 1388 and v\> ~ uq + v = 1412.5, we 
find v =24.5 cm _1 =3 meV. Also, the coupling translates into energy swinging between 
the 2 oscillators with frequency v : this is best seen by monitoring both angles after one 
of them has been changed slightly from its equilibrium value: they vary sinusoidally 
and out of phase with the expected periodicity, 1.4 ps. If Q\ and # 2 be the angles HCH, 
then, by symmetry, the amplitudes of the sym and a-s modes are simply £ = (#i + #2)/2 
and A = (0i — 6 2 ) /2. At low energy contents these quantities are perfect sinusoids with 
frequencies u a and z/j,, respectively; at higher contents, they are affected by amplitude 
modulation, whose frequency and depth obey (4) and (5). These relations will be used 
to study the changes in, and exchange between, the 2 modes with energy content. 

5. Case study: C2H4, energizing a single vibrational mode 

As an example, take as target mode the symmetric scissor mode. Normal mode analysis 
(provided by the semi-empirical codes) determines its wave number to be 1388 to the 
accuracy of this simulation. A periodic perturbation with T = 0.025 ps (1333 cm -1 ), 
making 5$ = —4%, was applied by setting the angles 8 between each pair of C-H bonds 
at 110 deg, from their equilibrium value, 114.6 deg: each such perturbation induces a 
potential energy increase of 0.9 kcal/mol or ~0.04 eV. Figure 4 shows the energizing 
curve: the excitation energy levels off near 110 kcal/mol, a limit imposed by the fact 
that the red-shifted proper frequency of the free motion has now fallen to 1335 cm" 1 , 
very near the excitation frequency so the energy gain per cycle tends to zero (see Sec. 
3). The peak angular excursion is also limited by stereoconstraints when the angles 
approach 180 deg: downward distorsion is indeed seen at the crest of the 9 (HCH) 
sinusoids. 

The excitation procedure was then repeated, starting from equilibrium but applying 
only 10 cycles to deliver 5 kcal/mol or about one mode quantum. The molecule was then 
left to execute its free motion in vacuum. No further evolution is visible for at least 100 
ps: on the screen, the molecule is observed to oscillate in an apparently pure symmetric 
scissor mode with constant amplitude. During this period, the relative energy content 
of the various modes is monitored in more detail by performing a FFT on the recorded 
total kinetic (or potential) energy. It is found indeed that the spectrum is dominated 
by a line very near the nominal target mode frequency. Its harmonics are more than a 
hundred fold weaker, indicating that the corresponding potential surface is very nearly 
parabolic. As expected in such circumstances, the mode is only very weakly coupled to 
other modes and none shows up in the spectrum above a relative level of 10~ 2 , except the 
i.p. symmetric CH bending (rock) mode at 835 cm" 1 and intensity e = 2 10" 2 relative 
to the scissor mode. A procedure to efficiently excite nearly single normal modes of 
vibration is thus demonstrated. 

At the same time, the coupling energies with other modes can be estimated using 
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Figure 4. Energizing curve for the symmetric scissor mode of C2H4. The periodic 
perturbation consists in imposing the same value, 110 deg, on both angles HCH, every 
0.025 ps . 



Sec. 2. The frequency separation of the bending mode from the target mode, 555 cm^ 1 , 
being so wide, one can assume that it is much larger than v in (5), so that the latter is 
approximately 10 cm -1 or 1.2 meV. This much is due to the geometrical similarity of 
the two modes. By contrast, the very close a-s scissor mode is not even detected because 
it holds a negative phase relation to the symmetric mode. Notwithstanding, it can be 
excited by increasing the energy content of the sym mode, but what happens then will 
best be understood after the following experiment is described. 

6. C2H4: anharmonicity, frequency dragging, mode locking and localization 

Here, the periodic perturbation (T = 0.025 ps) is not applied to both CH bond angles, 
as in Sec. 4, but to one of them only, No 1 (asymmetric excitation). As a result of their 
coupling, oscillator No 2 is also excited and they share the deposited energy, so that 
both scissor modes, E, A, are indeed excited (see beginning of Sec. 4). We followed 
their frequencies (fig. 5) as a function of the energy content of the molecule (varied 
by increasing the number of energizing cycles; see Sec. 3). It is observed that the 
anharmonic shift is initially small, then increases abruptly when the energy content 
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Figure 5. Illusrating mode freqency dragging, a and b represent the wave numbers of 
the symmetric and a-s modes, respectively, as a function of molecular energy content. 
Note the merging of the two curves at the higher energies, a phenomenon due to 
anharmonicity alone. 



nears 3 mode quanta. Ultimately, both frequencies merge when the energy exceeds 4 
quanta : this demonstrates the phenomenon of frequency dragging and mode locking (see 
Minorskyjl962 ). All the while, both amplitudes increase almost in the same proportion 
(although saturation is ultimately visible on the strongest mode). The behaviour of 9\ t 2 
illustrate still another molecular feature, namely localization (fig. 6): at low deposited 
energies, the dynamics is linear and both resonant oscillators share the energy equally 
but alternatively, according to Sec. 2; after an identical initial increase of both peak 
angles with input energy, the oscillator which was not initially excited starts to weaken 
and is ultimately completely quenched and the motion becomes confined to that angle 
which was initially perturbed. All these phenomena are continuous and linked with 
anharmonicity and saturation, and do not imply higher-order coupling terms (see end 
of Sec. 2). 
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Figure 6. 6\ and 02 represent the maximum angular excursions of the two HCH 

angles. Note the quenching of oscillator 2 by oscillator 1 at the higher energies. 



7. C 2 H 4 : locked "relaxation" oscillations 



Returning to the symmetric excitation of Sec. 5, we study the effects of increasing 
the energy content of the sym mode by the same procedure. There appears to be 
a threshold for the latter between 6 and 7 kcal/mol, beyond which, periodically, the 
a-s mode develops exponentially then dies out similarly. The turn occurs when the 
sym mode has provided so much of its energy that it falls below the threshold and 
therefore can no longer sustain the growth of the a-s mode (fig. 7). As behoves such 
relaxation mechanisms, the time constant of these evolutions decreases as the energy 
content increases. However, stereoconstraints (Sec. 5) progressively hinder the growth 
of the weaker (a-s) mode, so much as to finally quench it completely, as was the case in 
Sec. 6. 

Obviously, the coupling between these two modes is not of the "linear" type 
considered in Sec. 2. Numerical solutions on that basis, even including the first-order 
non-linear corrections, never exhibit the excitation of the a-s mode when, initially, only 
the sym mode is excited. That is, unless some sort of asymmetry is introduced, even 
though very weak, such as a slight inequality in z/ 12 or # 12 . One is therefore led to 
conclude that the growth of the a-s mode is made possible only by the unstable character 
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Figure 7. Mode locking and relaxation oscillations. With more or less delay following 
strong initial excitation of the sym scissor mode alone, the a-s mode is energized in 
bursts of decreasing amplitude as the energy content increases. The a-s mode wave 
number remains locked to that of the sym mode (see spectrum, fig. 8). Both behaviours 
are distinctly at variance with that of linear oscillators illustrated in fig. 1. At the 
limiting level of the energizing curve in fig. 4, the a-s mode is completely quenched. 



of the exciting mode at sufficiently high energy content. Such instabilities are known 
from non-linear mechanics to occur at the crossings of otherwise stable trajectories (see 
Wu[2005j). At these points, the trajectories are distorted so as to form two new separate 
trajectories, each connecting two different parts of the initial trajectories, and covering 
different energy domains. In molecular physics this is the well-known phenomenon of 
potential surface avoided crossing and repulsion, yielding a lower- and a higher-energy 
surface (see Herzberg[1945j). In the present experiment, we observe the swinging of 
the system between the two dips of the new lower-energy surface, which still essentially 
define the characters of the sym and a-s modes. The threshold energy for this to happen 
is obviously the energy at the tip of the intermediate hump between the dips. 

Beyond that, we expect both modes to have the same frequency, progressively red- 
shifted by anharmonicity as the energy content is increased. This is borne out by the 
spectrum of fig. 8, for instance, which corresponds to an energy content of 15 kcal/mol. 
The main peak of £ implies that this mode oscillates at 1388 cm" 1 and its 2 nearest 
side-bands, that it is modulated with frequency 8 cm -1 ( period 4.3 ps) as evidenced 



Molecular relaxation paths 



16 




v (cm" ) 



1370 1375 1380 1385 1390 1395 1400 1405 1410 



Figure 8. Spectra of the angular variation of £ and A, the sym (heavy line) and a-s 
(light line) scissor modes; see Sec. 8 for interpretation. 



by fig. 7. As to A, its 2 main peaks are symmetric with respect to the main peak of 
E, implying that this mode has the same frequency and the separation of the 2 peaks 
implies that it is also modulated, and with the same periodicity (cf. fig. 7). These 
frequency relations are maintained as the energy content varies: when the a-s mode is 
observable, it is locked to the symmetric mode. 

This behaviour cannot be emulated by coupled-mode theory unless a non-linear 
correction term in (x\ —X2) 2 is included. Such a term makes the system non-integrable, 
i.e. unstable. It is also responsible for the side-bands apparent in fig. 8. Indeed, the 
rate of decrease of these side-bands with distance from the center frequency is a measure 
of the time constant of the exponential rise and fall of the mode amplitudes. 



8. Higher-order resonances 

The next experiment is intended to illustrate "non-linear" (or "parametric" ) resonances, 
which give rise to combination frequencies and Fermi resonance. Here, the target mode 
is the sym CH stretch. An equal extension of 1.15 A (from the equilibrium value, 1.1 
A ) is applied periodically, every 0.011 ps (3000 cm -1 , 5 = 5%) to all 4 C-H bonds. 
After 10 cycles, the acquired excess energy is 9 kcal/mol (0.39 eV), nearly the energy of 
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Figure 9. Higher-order resonance. The symmetric CH stretch mode v\ is initially 
excited. Immediately thereafter, some of its energy begins to be transferred periodically 
to both the CC stretch, vi, and sym scissor, 2/3, modes, all 3 modes being of same 
symmetry species A g . Again, the envelope shapes witness to the non-linear character 
of the coupling giving rise to a 2nd-order resonance v\ = v-i + ^3. For about 10 ps, the 
spectra of these motions exhibit only these 3 frequencies. 



1 mode phonon. The free motion of the molecule is monitored thereafter. 

Almost immediately, the angles HCH and the bond length C-C are set in oscillation 
together with the C-H bond lengths (fig. 9). The former experience periodic surges in 
amplitude, for ~ 0.5 ps, every 2.1 ps, perfectly in tune with a strong decrease of the 
CH bond length amplitude of vibration. This suggests a very non- linear regime with 
energy swinging between the C-H motion, on the one hand, and the HCH and C-C 
motions on the other hand. The spectrum of each type of motion is indeed very pure 
, peaking respectively at 3239-3243, 1400-1403 and 1839-1840 cm" 1 . Obviously, the 
former is the sum of the latter two: this resonance occurs only because the molecular 
energy content is now high enough to perturb the parabolic shape of the potential energy 
surfaces to the point of inducing a significant higher term in the mode coupling terms 
(see end of Sec. 2). Because of the periodic energy exchange between modes, each band 
is accompanied by side-bands correspondig to sum and difference frequencies. Note 
that, due to imperfections of the AMI type of chemical modeling used here, the wave 
numbers differ slightly from experimental ones, without affecting the phenomenological 
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Figure 10. Follow-up to fig. 9. About 10 ps after energizing the CH stretch mode, 
the spectra become more complicated. The spectrum of the total molecular kinetic 
energy represents here the new (and final) energy distribution between modes, in 
kcal/mol/cm -1 . The initial sym CH str energy is now partitioned between this mode 
and the a-s CH str mode (58 %), the CC str mode (21 %) and the 2 scissor modes (20 
%), totalling 99 % between them. 



behaviour. 

Further analysis of the bond angle and length variations reveals that the first band 
corresponds to pure symmetric C-H stretch (all 4 bonds being sychronized) and the 
second to pure symmetric scissor, to the accuracy of the present simulation. 

After about 10 ps, the motion becomes more complex: fig. 10 shows the FFT of the 
total molecular kinetic energy. In more detail, the 2 pairs of C-H stretches are no longer 
synchronized, neither are the 2 HCH angles. In order to break their motions down into 
modes, we take the spectra of (CiH 3 +CiH4)±(C 2 II 5 +C2H 6 ) and H 3 CiH4±H 5 C2lI 6 . We 
now have 2 bands for each of these combinations: respectively 

sym CH str 3210, 3245 cm" 1 , 

a-s CH str 3220, 3264 cm" 1 , 

sym scissor 1377, 1413 cm -1 , 

a-s scissor 1389, 1432 cm -1 . 

The CC str is only slightly shifted to 1832 cm -1 . Each of the 4 CH str wave 
numbers is found to be the sum of the wave numbers of the CC str and one of the two 
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Figure 11. Illustrating irreversible transition from a regular motion to a partially 
chaotic regime, following the stretching of a C-H bond in ethylene (see Sec. 9) 



scissors. Thus, the same type of non-linear resonance now extends to the scissor and 
CH str modes (of course, the 2 other CH str modes can hardly couple with the scissors 
to generate a resonance). 

Here again, the qualitative aspects of these observations can only be mimicked, 
using the equations of Sec. 2, if the term {x\ — x 2 ) 2 is included, which allows the formation 
of combination frequencies. The other term, affecting explicitely the oscillators' resonant 
frequencies, brings no qualitative change to the picture. A quantitative determination, 
by cut-and-try, of the perturbative terms responsible for specific observations is outside 
the scope of the present work 

Finally, after about 100 ps, steady-state seems to be established, and the relative 
amounts of energy in the scissor, CC str and CH str types of motion are 0.34, 0.36 and 1, 
respectively; the rest is negligible. Thus, only 5 modes out of 12 have been excited, and 
then, not even equally. A similar conclusion can be drawn from the selective excitations 
of other modes. For instance, if one mode quantum is deposited in the CC stretching 
vibration, less than 1% is transferred to the other modes. 
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9. Transition to irreversible IVR 

When the energy deposited locally or in a single mode becomes large, new phenomena 
are observed. In the following experiment, one of the C-H bonds was stretched to 1.2 
A , giving the molecule an excess energy of 20 kcal/mol and left to move freely. Figure 
11 shows the autocorrelation of the total kinetic energy as a function of time. The 
first picoseconds are characteristic of regular (quasi- harmonic) motion: note the strong 
periodic peaks with a weak and slower amplitude modulation superimposed, indicating 
the dominance of CH str oscillations. Soon, a steep fall in the peak amplitude becomes 
apparent and the pattern changes condiderably: the correlation is substantially reduced 
and it becomes difficult to recognize any regularity in the modulation of its envelope. 
At the same time, many modes other than the stretches show up in the spectra. This 
experiment illustrates an irreversible path to chaos, namely increasing energy density. 
However, the deposited energy here is soon diluted, even in so small a molecule, and 
the irreversible process comes to a stop. 

10. Large molecules 

Another favorable factor for chaos is the size of the oscillating system, because this 
increases the spectral density of modes, and thus favours the overlap of resonances 
(Ford [1992], Sibert et al [1984]) and opposes the spatial dilution of energy. In order 
to illustrate this phenomenon, the same type of computational experiments were 
performed on a much larger hydrocarbon molecule, composed of interlinked peri- and 
catacondensed carbon rings, and including benzene, pyrene and pentagonal rings, 8 
oxygen, 2 sulphur and 1 nitrogen atoms, 11 tertiary, 44 secondary and 9 primary H 
atoms for a total of 101 atoms, fully described in Papoular |2002j . The structure is 3-D 
and has linear extensions of order 1 nm. The lengths and relative orientations of bonds 
are initially computationally optimized, at rest, for a minimum total potential energy. 

The large number of normal modes in this case (297) uncovers new interesting 
behaviour, which extends over longer periods of time. The computational simulation 
then is no longer feasible at the higher-quality level, so a less sophisticated software 
was used here: Molecular Mechanics MM2. While not affecting severely the general 
features of the underlying physics, this is very much faster and allowed the extension of 
simulations up to 1 ns (Papoular 2002J). 

Another difficulty lies in devising a priori a system of simultaneous perturbations of 
all bonds so as to emulate a prescribed normal mode. Fortunately, this turns out not to 
be essential for present purposes. In the following experiment, the target modes are the 
CH stretches. The energizing procedure consists in imposing periodically ( T = 0.012 
ps, v = 2800 cm" 1 ) to 2 CH bonds at opposite edges of the molecule (one aliphatic and 
one aromatic) an extension of 1.15 A , from their equilibrium value, ~1 A . Twenty 
cycles are enough to deposit 8.5 kcal/mol (~1 mode quantum) in the molecule. 

The following free-motion relaxation is not smooth. Its different phases are 
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Figure 12. Relaxation of a 101-atom hydrocarbon molecule following selective 
excitation of CH sretchings, illustrated by the evolution of the standard deviation 
of the total kinetic energy (kcal/mol) . Note the abrupt downward transitions near 50 
and 170 ps. 



strikingly illustrated by monitoring the standard deviation, o~k, of the total kinetic 
energy, E K : statistical arguments show that, roughly, o~k/E k = n -1 / 2 , where n is 
an effective number of vibrational modes, considered as gas particles (see Reif[1965j). 
Figure 12 shows ok (computed over 8 ps for each point) as a function of time. This 
uncovers 3 different phases separated by "abrupt" transitions and differing in the 
spectral distribution of energy: 

a) to 50 ps. The energy remains localized mainly in the two perturbed bonds; 
the spectrum of Ek exhibits only their two different local frequencies, 2775 and 2809 
cm -1 , still very close to the excitation frequency. This spectral sparseness is the cause 
of the large fluctuations in o~ K . 

b) 50 to 160 ps. ok now decreases smoothly because many other C-H bonds have 
been strongly excited. They are mainly aliphatic and oscillate at ~2820 cm" 1 , redward 
of their natural frequency because of their high energy content. Some energy has gone 
into in-plane CH bendings and CC stretches. 

c) 160 to 1000 ps. In this final phase, derealization is complete, and both aliphatic 
and aromatic CH str are found near their natural frequencies, 2835 and 2940 cm" 1 
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Figure 13. Average number of mode quanta, his, per mode v, (solid line), deduced 
from FFT of the total kinetic energy over the final 8 ps of the dynamics, for the case 
of fig. 12. The dashed line represents the thermodynamic equilibrium at temperature 
130 cm -1 , computed for the same energy content, 2900 cm -1 . 



(to the accuracy of this software). Most importantly, lower frequency modes are 
progressively and smoothly excited down to the lowest. Figure 13 shows the spectral 
energy distribution during the last 8 ps, in terms of average number of mode quanta per 
mode. 

We take this to be steady-state distribution of energy because fig. 12 shows that 
ok is close to its asyptotic value, estimated at 0.36. This value implies that ~140, out of 
the total 297, modes have been effectively excited, most of them at the lower end of the 
spectrum, mainly skeleton modes. Of the energy initially deposited in the molecule, 28 
% are still in the target CH str modes. It must be stressed that these are average values 
for a steady but dynamic state. Energy is constantly exchanged between modes, and this 
corresponds, in quantum mechanics, to the flow of energy throughout all the molecular 
states of the portion of phase space that is available at this energy content, taking into 
account selection rules and all other restrictions. If the 297 quantum coupled oscillators 
had somehow settled down into a statistical thermodynamic equilibrium with the same 
energy content, its temperature would have reached 130 cm -1 (see Allamandola et al 
[1989J), and the average number of mode quanta in mode v would be l/[exp(u/130) — 1], 
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Figure 14. Same as fig. 12, following deposition of 810 cm -1 in a CC str mode. Here, 
the vertical scale is logarithmic so as to highlight the initial exponantial character of 
the relaxation. 



as illustrated by the dashed curve in fig. 13. The difference is striking and due to 
the much larger phase space covered in a fully ergodic motion, subject only to energy 
conservation, which greatly favours low frequencies and penalizes high frequencies. 

The total relaxation time is about 1 ns. By contrast, the transition time is only 
of order 10 ps. During this time, larger and larger parts of the molecule are seen to be 
quickly set in motion. This signals again the existence, at high local energy densities, of 
metastable states which soon lead to significant redistribution of energy among modes 
and changes in the spectra. 

Different types of relaxation are encountered when initially energizing different 
normal modes. In the following and last experiment, one of the CC bridges between 
aromatic ring clusters was energized by periodic perturbations of its length, with 
frequency 926 cm -1 . Again, the frequency selectivity of the process makes for a very 
sparse spectrum in the first 10 ps of the free motion: only 2 modes are excited, at 876 and 
948 cm -1 , which oscillate coherently (i.e. exchange energy periodically, as in Sec. 2), 
giving rise to beatings in the total kinetic and potential energies, with period 0.225 ps, 
corresponding precisely to their frequency difference. However, the relaxation process 
immediately sets in. It is now exponential for about 200 ps, with an exponent of 135 
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Figure 15. Same as fig. 13, for the case of fig. 14. The dashed line is for 
thermodynamic equilibrium at T = 70 cm -1 , corresponding to the same energy 
content, 810 cm -1 . 



ps (fig. 14). Afterwards, it begins to level off but still exhibits long period fluctuations, 
also seen in the auto-correlation functions of the total energies. This indicates bulk 
structural vibrations at very low frequency, as borne out by the final spectral energy 
distribution, fig. 15, where the dominance of low frequencies is even more strongly 
pronounced than in the previous case. 

11. Conclusion 



Single-mode (or nearly so) excitation of isolated molecules allows one to explore 
different relaxation paths towards different dynamic steady-states. Along the way, it 
is possible to indentify bottlenecks and determine coupling strengths between modes. 
The high temporal resolution of computational modeling uncovers phenomena which 
are characteristic of the non-linear domain of molecular dynamics, but not very easy to 
demonstrate experimentally, such as local modes, frequency dragging and mode locking, 
as well as metastable states (saddle points) and relaxation oscillations. 

Coupled-oscillators theory helps understanding the underlying physics, disentan- 
gling complicated states and quantifying the parameters responsible for these phenom- 
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ena. 

The final relaxed states of the molecules are found to be different from each other 
and from thermodynamic equilibrium. Computational experiments may therefore help 
preparing molecules for specific chemical reactions, or predicting their IR emission 
following specified excitation, which was the initial goal of this research. 
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